library(xCell)
library(tidyverse)
library(ggpubr)

# gtex database --------
gtex <- read_delim('mission/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct.gz',
                   skip = 2)
gtxcell <- gtex |>
  distinct(Description, .keep_all = TRUE) |>
  select(-1) |>
  column_to_rownames('Description') |>
  as.matrix() |>
  xCellAnalysis()

gtxcell |>
  as_tibble(rownames = 'celltype') |>
  filter(celltype == 'Plasma cells') |>
  pivot_longer(-1) |>
  filter(!str_detect(name, 'EBV')) |>
  slice_max(value, n = 15) |>
  ggplot(aes(fct_reorder(name, value), value)) +
  geom_col() +
  coord_flip() +
  theme_classic(base_size = 12) +
  labs(title = 'Enrichment of plasma cell in healthy tissue\n by xCell',
        subtitle = 'GTEx data',
        y = 'Enrichment score',
        x = 'Tissue')

gtxsig <- gtxcell |>
  as.data.frame() |>
  head() |>
  select(Pancreas) |>
  as.matrix() |>
  xCellSignifcanceBetaDist()

gtxsig |> head()

# HPA database ---------
hpca <- read_delim('mission/rna_tissue_hpa.tsv')

unique_ensem <- hpca |>
  distinct(`Gene name`, .keep_all = TRUE) |>
  pull(Gene) 

hpaxcell <- hpca |>
  filter(Gene %in% unique_ensem) |>
  pivot_wider(id_cols = `Gene name`, names_from = 'Tissue', values_from = 'nTPM') |>
  column_to_rownames("Gene name") |>
  xCellAnalysis()

hpaxcell |>
  as_tibble(rownames = 'celltype') |>
  filter(celltype == 'Plasma cells') |>
  pivot_longer(-1) |>
  slice_max(value, n = 15) |>
  ggplot(aes(fct_reorder(name, value), value)) +
  geom_col() +
  coord_flip() +
  theme_classic(base_size = 12) +
  labs(title = 'Enrichment of plasma cell in healthy tissue\n by xCell',
       subtitle = 'HPA data',
       y = 'Enrichment score',
       x = 'Tissue')

hpaxcell |>
  as_tibble(rownames = 'celltype') |>
  filter(celltype == 'Plasma cells') |>
  pivot_longer(-1) |>
  filter(str_detect(name, 'marrow'))

xCell.data$signatures |>
  as.list()

plasma_index <- xCell.data$signatures |>
  names() |>
  str_subset('lasma')

pc_xcell <- xCell.data$signatures[plasma_index]

pc_gene <- pc_xcell |>
  map(\(x)x@geneIds)

pc_kegg <- pc_gene |>
  as.list() |>
  reduce(union) |>
  WebGestaltR(enrichDatabase="pathway_KEGG", interestGene = _,
              referenceSet = 'genome',
              interestGeneType="genesymbol", isOutput=FALSE)

# Tabula sapiens scRNA----------
source('00_util_scripts/mod_bplot.R')
  
tabula <- read_csv('mission/Tabula_Sapiens_metadata.csv')

tabula |> count(cell_ontology_class, organ_tissue) |>
  group_by(organ_tissue) |>
  calc_frac_conf_on_grouped_count() |>
  filter(str_detect(cell_ontology_class, 'plasma cell')) |>
  ggplot(aes(fct_reorder(organ_tissue, fraction), fraction, fill = fraction)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'none') +
  labs(title = 'Fraction of plasma cell in healthy tissue',
       subtitle = 'Tabula scRNA-seq',
       y = 'Fraction of cells',
       x = 'Tissue')
